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We present a new adaptive resolution technique for efficient particle-based multiscale molecular 
dynamics (MD) simulations. The presented approach is tailor-made for molecular systems where 
atomistic resolution is required only in spatially localized domains whereas a lower mesoscopic level 
of detail is sufficient for the rest of the system. Our method allows an on-the-fly interchange between 
a given molecule's atomic and coarse-grained level of description, enabling us to reach large length 
and time scales while spatially retaining atomistic details of the system. The new approach is 
tested on a model system of a liquid of tetrahedral molecules. The simulation box is divided into 
two regions: one containing only atomistically resolved tetrahedral molecules, the other containing 
only one particle coarse-grained spherical molecules. The molecules can freely move between the two 
regions while changing their level of resolution accordingly. The coarse-grained and the atomistically 
resolved systems have the same statistical properties at the same physical conditions. 



I. INTRODUCTION 



Many problems in complex soft matter systems 
are inherently multiscale in nature, i.e., micro- 
scopic interactions are coupled strongly to meso- 
and macroscopic properties. Despite the increas- 
ing computational power and ongoing efforts to en- 
hance the efficiency of molecular dynamics (MD) 
integration algorithms^, EI [J [3 all-atom 
MD simulations are often incapable to cover the 
time and length scales needed in order to reach re- 
laxation in a typical molecular system, such as a 
polymer solution or melt. In many cases it is also 
questionable, whether the huge amount of detail 
information might not even obscure the relevant 
structural information. On the other hand, details 
of the chemistry do not affect universal power laws 
but the prefactors of these power laws, which can 
vary by several orders of magnitude themselves. 
Thus even on the more coarse-grained level it is 
advisable to keep a link to the specific chemistry 
under investigation^, Q- In addition typical soft 
matter systems can be quite inhomogeneous in a 
way, that different regions within one system are 
sufficiently described by more or less detail. A 
consistent and at the same time highly efficient 
ansatz to understand modern soft matter systems 
(both of synthetic as well as biological origin) has 
to take such considerations into account. One first 
way to tackle this, was to reduce the number of de- 
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grees of freedom by a systematic coarse-graining, 
which retains only those degrees of freedom that 
are relevant for the particular property of interest. 
Examples of molecular systems where the coarse- 
graining approach has been used with much suc- 
cess are fluids 0, lipid biWersJflJll El El 13 , 
and polymer systems [TR lifil Il7t llSj . 

Since some specific chemical details are usually 
lost in the coarse-graining procedure, much effort 
has been devoted recently to the development of 
multiscale modeling approaches, where different 
parts of the system are modeled at different lev- 
els of detail to account for the local resolution 
requirement [H HJ M, El H . In the dual- 
resolution modeling approach for stud ying the be- 
havior of polymers near metal surfaces p3 E^ . l24j . 
for example, polymer chain ends that interact with 
the metal surface are represented partially atom- 
istically while the remaining parts of the polymers, 
where lower resolution is adequate, are represented 
as bead-spring chains. However, the switching be- 
tween different levels of resolution, i.e, the atomic 
and mesoscopic, is not allowed during the course of 
that MD simulation, and therefore the initial level 
of detail and thereby the number of degrees of free- 
dom in the system remain unchanged. Since the 
chain ends that interact with the metal typically 
remain close to the surface and they only contain 
a small fraction of the whole system, an adaptive 
on-the-fly change of the molecules' resolution is 
not strictly required for this class of systems. An- 
other approach reported in the literature concerns 
the link between quantum mechanical and classi- 
cal MD simulations. In this QM/MM approach 
a small subset of the system is defined and con- 
sidered as quantum mechanical, while the rest is 
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treated by a classical force field simulation. Here 
the atoms as well as the regions of the different 
regimes are fixed from the very beginning |25j| , re- 
stricting the application to rather specific cases. 

In contrast, MD simulations, in which the spa- 
tially localized atomistic domains frequently ex- 
change particles with the remaining mesoscopic 
part of the system, would allow for much wider 
applications. Then it would be possible to define 
certain areas or develop criteria for certain situa- 
tions, which ask for a more detailed view, while 
the rest of the system can be treated on the more 
coarse-grained level. It is the purpose of this pa- 
per to present a first attempt of an MD simula- 
tion of that kind. Because of that we here at first 
restrict ourselves to a most elementary model sys- 
tem. Of course one also could resort to Monte 
Carlo simulations, where different levels of detail 
are combined. This actually would be somewhat 
simpler because of the purely stochastic nature 
of this simulation method. Since we are eventu- 
ally aiming at molecular systems, where collec- 
tive motions are crucial, we decided to stick to 
the MD approach. Existing hybrid MD meth- 
ods that concurrently couple different length scales 
have been developed to study solid state systems, 
where atomistic MD was either combined with 
the finite elements methodpCjl I2H I2II or it was 
linked to a quantum mechanical model|29(. To our 
knowledge, however, in pursuit of this objective 
no adaptive hybrid atomistic/mesoscale particle- 
based MD method, which would allow to dynam- 
ically adjust the level of detail, which means the 
adjustment of the degrees of freedom in the sys- 
tem, has been developed so far. 

In this paper we present a novel adaptive reso- 
lution MD scheme that combines a full atomistic 
description of a desired region of the system with 
a mesoscale treatment of the remaining part. The 
key feature of the new method is that it allows to 
dynamically adapt the level of a given molecule's 
resolution according to its position in the system. 
Hence, the number of degrees of freedom is not a 
conserved quantity in our MD simulations. Fur- 
thermore, the presented method is not restricted 
to couple only the atomistic and mesoscopic levels 
of detail but can also be applied to systems with 
mesoscopic domains that are described at differ- 
ent levels of coarse-graining. Here we present a 
first test case, showing that such an approach is 
feasible. Therefore the new approach is tested 
for a simple model liquid consisting of tetrahe- 
dral molecules. The simulation box is divided 
into two regions: one containing " atomistically" 



resolved tetrahedral molecules, the other contain- 
ing coarse-grained spherical molecules. Molecules 
are allowed to freely move between the two regions 
while changing their level of resolution accordingly. 
The results show that the statistical properties of 
the corresponding fully atomistic system are accu- 
rately reproduced by using the proposed hybrid 
scheme. In particular, gradients in the chemi- 
cal potential across the artificial interface where 
the resolution changes and corresponding spurious 
fluxes can be avoided. Although we applied the 
new method here to a generic test system it should 
find application in more realistic physical systems. 
There has been an initial attempt reported in Ref. 
|3(ij | , where MD simulation of an inhomogeneously 
coarse-grained system of liquid methane has been 
described. There the starting point were two es- 
tablished models for methane, a five site atomistic 
and a one site spherical where the interaction be- 
tween molecules of different species was derived by 
standard Lennard- Jones mixing rules with the hy- 
drogens of the atomistic invisible to the spherical 
molecules. In this case it turned out that an ef- 
fective flux between the two different regimes oc- 
curred. This effect is partially due to the differ- 
ent equilibrium state points described by the two 
models, but in any case in that approach no ef- 
fective potential between coarse-grained molecules 
was derived. Our approach differs from that since 
here we derive such potentials in such a way that 
the two regimes are in true thermodynamic equi- 
librium. 

The organization of the article is as follows. In 
section 2 the methodology is presented, whereby 
in the first step a coarse-grained model is derived 
and parameterized from a fully atomistic system 
and then, in the second step, the atomic and meso- 
scopic length scales are systematically coupled in a 
hybrid atomistic/mesoscopic model. The compu- 
tational details are given in section 3. The results 
and discussion are presented in section 4, followed 
by conclusions in section 5. 



II. METHODOLOGY 



In this section, the model systems are described 
and the methodology for the adaptive multiscale 
MD simulations is presented. 
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Models 

All-Atom Model 

First, we introduce our reference explicit all- 
atom (ex) model. Consider a system of n tetra- 
hedral molecules consisting of N = 4 atoms of the 
same mass mo connected by anharmonic bonds as 
presented in figure ^ (a) (consider only the right 
red molecule). 




(b) 

FIG. 1: (a) The on-the-fly interchange between the 
atomic and coarse-grained levels of description. The 
middle hybrid molecule is a linear combination of 
fully atomistic tetrahedral molecule with an addi- 
tional center of mass particle representing the coarse- 
grained molecule, (b) Snapshot of the hybrid atom- 
istic/mesoscopic model at p* — 0.1 and T* = 1 
(LJ units). The red molecules are the explicit 
atomistically resolved tetrahedral molecules, the blue 
molecules are the corresponding one particle coarse- 
grained molecules. 



All atoms in the system interact according to 
a purely repulsive shifted 12-6 Lennard- Jones po- 
tential with a cutoff at 2 1 / 6 <r: 

/M(^) 12 -fe) 6 + !]; ^<2V ff 

\ 0; r lajP > 2V6 CT 

(1) 

where ri a jp is the distance between the atom ia of 
the molecule a and the atom j/3 of the molecule (3. 
We use e as a unit of energy. All atoms have the 
same excluded volume diameter a, where a is the 
unit of length. The neighboring atoms in a given 
molecule a are linked via an attractive FENE po- 
tential 

jjatom ( \ _ 

| -ifcR§ln[l - (^ff ) 2 ]; r iaja < R ^ 

>Ro 

with divergence length Ro = 1.5a and stiffness 
k = 30e/a 2 , so that the average bond length is 
approximately l.Oer for fc^T = e, where T is the 
temperature of the system and ks is Boltzmann's 
constant. The functional form of these potentials 
and their parameters are the same as usually em- 
ployed in polymer MD simulations|l5|. 

Coarse-Grained Model 

Next, we map the atomistic model to a coarse- 
grained (eg) mesoscopic model. For the latter we 
have chosen a system composed of n one-particle 
molecules schematically depicted in figure ^ ( a ) 
(consider only the left blue molecule). A given 
coarse-grained molecule a in the system has a mass 
M a = 4mo equal to the total mass of the explicit 
tetrahedral molecule. All rotational and vibra- 
tional degrees of freedom of atomistically resolved 
tetrahedral molecules are thus removed, and the 
number of nonbonded interactions is strongly de- 
creased as well. 

We shall now find an effective pair potential 
between coarse-grained molecules such that the 
structural properties of the underlying atomistic 
model are reproduced. There is usually no unique 
way to coarse-grain to an effective pair potential, 
which is in general temperature and density depen- 
dent |3lL l32l l33 | . just as in statistical mechanics 
there are different ways to perform a renormaliza- 
tion group step. Here, we extract the effective pair 
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potential from a center-of-mass radial distribution 
function (RDF cm ) of the reference atomistic model 
using the potential of mean force PMF(r) as 

U cm (r) « PMF{r) = -k B T log 5e c ™(r), (3) 

where 3e™( r ) i s the RDF cm of the all-atom sys- 
tem and U cm (r) is the derived effective pair 
potential j34[. The effective potential obtained in 
this way is correct only in the limit of zero den- 
sity, where the many-body contributions vanish. 
For systems with nonzero densities in principle 
many-body interactions would be needed just as 
for the classical renormalization group theory in 
statistical mechanics. In a similar spirit we here 
resort to the (expected) relevant part in order to 
gain a significant speed up in our simulations. Be- 
cause of that we use the PMF as the initial guess 
for the effective pair potential in systems with 
nonzero density. Then this is further fine-tuned 
until the RDF cm s and pressures of the reference 
atomistic and coarse-grained systems match [34) . 
As it turns out (see the Results and Discussion 
section) the effective pair potential acting between 
our coarse-grained molecules is significantly softer 
than the pair potential between atoms of the re- 
solved molecules in accordance with the results 
previously found in the literature |33l l35| . 

Transition Regime: Hybrid Atomistic /Mesoscopic 
Model 

Let us now introduce a hybrid explicit/coarse- 
grained (ex-cg) model. Consider a box of n 
molecules where one half of the box is occupied 
by atomistically resolved four-atom tetrahedral 
molecules while the other half is filled up with the 
same number of corresponding coarse-grained one- 
particle molecules as schematically presented in 
figuren(b). The two domains exchange molecules 
which adapt their level of resolution accordingly. 
To ensure that the transition between the two dif- 
ferent levels of description is smooth, i.e., the ro- 
tational and vibrational degrees of freedom of a 
given molecule are gradually 'switched on/off' as 
it crosses the boundary between the atomistic and 
coarse-grained domains, we also introduce an in- 
terface layer for handshaking between atomistic 
and mesoscale regions. In this regime also the ro- 
tational and vibrational velocities have to be rein- 
troduced in a consistent way. Due to the periodic 
boundary conditions employed in our simulations 
there are in fact two such layers as depicted in 
figure 03 (b). 



Since the total number of molecules n in the sys- 
tem is a conserved quantity in our constant tem- 
perature simulations, we sample the phase space 
from the nVT ensemble. However, the total num- 
ber of degrees of freedom is not constant in this 
model. 

An alternative way compared to the similari- 
ties with renormalization group methods, which 
probably describes the situation even better, is the 
comparison to a first order phase transition. The 
rotational and vibrational part of the free energy 
then can be viewed as the latent heat at this transi- 
tion. At equilibrium, conditions analogous to two 
phase coexistence, 

f-^ex l^cgi Pex — Peg: T ex — Tcgi (4) 

must be automatically satisfied, where fi ex , p ex , 
T ex and (j, cg , p cg , T cg are the chemical poten- 
tials, pressures, and temperatures of the liquid in 
the atomistic and coarse-grained domains, respec- 
tively. These conditions (JIJ assure that there is no 
net flux of molecules between the atomistic and 
coarse-grained regions. To keep this absolute re- 
quirement also then defines a central task. This 
guarantees that the liquid is homogeneous across 
the box as it is in the reference fully atomistic sys- 
tem. From a molecular point of view, the artificial 
resolution boundary must be essentially invisible, 
i.e., the molecules have to cross the border with- 
out experiencing any 'barrier'. Our approach to 
reach this objective is presented in the proceeding 
subsection. 



Adaptive Resolution Scheme 

To allow a coarse-grained molecule to find 
an energetically permissible orientation with re- 
spect to its neighboring molecules (when it leaves 
the coarse-grained domain and is remapped into 
the atomistically resolved four-atom tetrahedral 
molecule) we introduce an interface layer between 
the atomistic and coarse-grained regions, which 
contains 'hybrid' molecules as presented in figure 
^ Each hybrid molecule schematically shown in 
figure ^a) (consider the middle molecule) is com- 
posed of a tetrahedral molecule with an additional 
massless center-of-mass particle serving as an in- 
teraction site. This is similar to the flexible TIP4P 
water model[36j where apart from the interaction 
sites on the three atoms of a water molecule an 
additional interaction site is introduced along the 
symmetry axis between the hydrogen and oxygen 
atoms. 
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Thus, each time a coarse-grained molecule a 
leaves that domain and enters the interface layer, 
it is remapped first into a hybrid molecule with 
the same center-of-mass position and a random 
orientation in space, where the relative positions 
of the tetrahedral atoms are taken from a molecu- 
lar configuration corresponding to a randomly cho- 
sen molecule from the atomistic regime. Each of 
the four explicit tetrahedral atoms in the hybrid 
molecule gains at this remapping a velocity equal 
to the velocity of the corresponding coarse-grained 
molecule to maintain the linear momentum of the 
molecule. In addition, the tetrahedral atoms are 
also assigned rotational/vibrational velocities cor- 
responding to atoms of a random molecule from 
the atomistic region, where we subtract the to- 
tal linear momentum of the latter molecule. In 
this way we ensure that the kinetic energy is dis- 
tributed among all degrees of freedom according 
to the equipartition principle as fceT/2 of aver- 
age kinetic energy per quadratic degree of freedom 
while retaining the linear momentum of the coarse- 
grained molecule. The center-of-mass interaction 
site moves obeying the constraints: 



R Q — 

V Q = 



Zjjq m ia v ic 

M a 



(5) 
(6) 



where R Q is a center of mass of the molecule a, 
Ti a is the position vector of the explicit tetrahe- 
dral atom ia in the molecule a, V a is the center- 
of-mass velocity of the molecule a, Vi a is the ve- 
locity of the explicit tetrahedral atom ia, and 
M a = J2ia m ia i s t ne total mass of the molecule 
a. In our case, mi a — mo and M a = 4mo for 
all i a = 1, . . . ,4 and a — 1, . . . , n. Each time a 
hybrid molecule crosses the boundary into atom- 
istic regime it is remapped into a four-particle 
tetrahedral molecule with the four tetrahedral 
atoms retaining their current velocities and po- 
sitions. In this model also the explicit tetrahe- 
dral molecules have the center-of-mass interaction 
sites, but only for the interactions with the hybrid 
and coarse-grained molecules. Of course, deep in 
the atomistic region, where the atomistically re- 
solved molecules do not interact anymore with the 
hybrid and coarse-grained molecules the center-of- 
mass interaction site can be omitted. Every time a 
four-particle tetrahedral molecule leaves the atom- 
istic region and enters into the transition regime 
it is mapped into a hybrid molecule with the four 
tetrahedral atoms retaining instantaneous veloci- 
ties and positions with the center-of-mass interac- 



tion site moving according to Eqs. (JjjJ and JSJ. 

Similarly, as a hybrid molecule crosses a bound- 
ary to the coarse-grained region it is mapped into 
a coarse-grained molecule with a velocity equal to 
the center-of-mass velocity of the hybrid molecule 
given by Eq. 

To couple the atomic and mesoscopic length 
scales we define in the spirit of thermodynamic 
perturbation approach|37l l38l ] the total inter- 
molecular force acting between centers of mass of 
molecules a and (3 as 



* a/3 — 

w{X a )w{X p )F a a t p + [1 - w(X a )w(Xp)]F<^, 



where 



Fatom 
a(3 



* iaj/3 

ia,j/3 



(7) 



(8) 



is the sum of all pair atom interactions between 
explicit tetrahedral atoms of the molecule a and 
explicit tetrahedral atoms of the molecule (3 and 



-patom 



Fcm 
a/3 



dU° 



Qjjcm 



<9R, 



a/3 



(9) 
(10) 



The vector Yi a jp = Ti a — Tjp is the relative posi- 
tion vector of atoms ia and j/3, R a /3 = R Q — R^ is 
the relative position vector of the centers of mass 
of the molecules a and [3, X a and Xp are the x 
center-of-mass coordinates of the molecules a and 
/?, respectively, and w is the weighting function 
that determines the 'identity' of a given molecule. 
The weighting function w G [0, 1] is defined in such 
a way that values < w < 1 correspond to a hy- 
brid molecule with extreme cases w — 1 and w = 
corresponding to a four-atom tetrahedral molecule 
and one-particle coarse-grained molecule, respec- 
tively. Hence, as soon as one of the two interact- 
ing molecules a and (3 is a coarse-grained molecule 
with no explicit tetrahedral atoms w(X a )w(Xp) = 
and F Q/3 = F^. 

We propose the following functional form of the 
weighting function w: 



w(x) 



cos' 



sin 



1 



d)\ 
d)} 
d)} 



d<x<%-d 
-f + d < x < -d 
—d < x < d 

f ~d<x< f 
-I < x < -# +d 



(11) 



6 



where a is the box length and d the half-width 
of the interface layer. The weighting function w 
is shown in figure [21 Our choice, which takes into 
account the periodic boundary conditions, is a par- 
ticularly simple way to ensure an interpolation be- 
tween w = and w = 1 that is monotonic, con- 
tinuous, differentiable and has zero slope at the 
boundaries to the atomistic and coarse-grained re- 
gions. Apart from these requirements we consider 
the precise functional form as immaterial. 



s 0.5 




o 1 — * ' 1 ^ — 1 1 1 ' 1 1 

-0.5 -0.4 -0.3 -0.2 -0.1 0.1 0.2 0.3 0.4 0.5 

x/a 

FIG. 2: The weighting function w(x) £ [0, 1] defined 
by Eq. ill U . The values w = 1 and w = correspond 
to the atomistic and coarse-grained regions of the hy- 
brid atomistic/mesoscopic system with the box length 
a, respectively, while the values < w < 1 correspond 
to the interface layer. Shown is the example where the 
half-width d of the interface layer is d = a/10. The 
vertical lines denote the boundaries of the interface 
layers. 

Exploiting the analogy with the quantum me- 
chanical mixed state description (3^, one can con- 
sider a hybrid molecule in the interface layer 
as a normalized linear combination of a four- 
atom tetrahedral molecule and a corresponding 
one-particle coarse-grained molecule. As a given 
molecule moves from the coarse-grained boundary 
of the interface layer to the atomistic boundary, 
w is gradually changed from to 1 and a coarse- 
grained molecule with only 3 translational degrees 
of freedom gradually turns into an atomistically 
resolved molecule with additional 3N — 3 = 9 ro- 
tational and vibrational degrees of freedom and 
a defined spatial orientation. The continuous, 
i.e., not instantaneous, 'identity' transition is re- 
quired since a hybrid molecule is given a random 
orientation at the coarse-grained boundary, and 
there can consequently be overlaps of its tetra- 
hedral atoms with the atoms of the neighboring 
molecules. Since at this boundary w = and 



the repulsive potential JJ atom given in Eq. Q is 
capped (see the Computational Details section), 
the forces acting on these atoms cannot diverge, 
however. While moving towards the atomistic re- 
gion the hybrid molecule slowly adapts its orienta- 
tion via the gradually increasing atomistic interac- 
tions with the neighboring molecules. Likewise, as 
presented in figure^(a), as a given molecule moves 
from the atomistic boundary of the interface layer 
to the coarse-grained boundary, w is continously 
changing from 1 to 0, the fully atomistically re- 
solved molecule gradually turns into the one par- 
ticle coarse-grained molecule while omitting all ro- 
tational and vibrational degrees of freedom and 
orientation. 

To elucidate the definition of the force calcula- 
tion in our model we rewrite Eq. J7J as 

F a /3 = 

w(X a )w(X p )F^ + [1 - w(X a )][l - w(X p )]F c ™ 

+ [l-w{X a )]w{X )¥^+w{X a )[l-w{Xp)]¥^. 

(12) 

From Eqs. JSJ and (|12l) one can then deduce that 
the pair force in Eq. J7J) is defined in such a way 
that two given atoms ia and jj3 in given explicit 
molecules a and (3 (w(X a ) = 1 and w(X/3) = 1) 
interact via the atomistic potential defined by Eqs. 
Q and J2J while two coarse-grained molecules 
(w(X a ) — and w(Xp) = 0) interact via the cor- 
responding effective pair potential U cm . Further- 
more, the coarse-grained molecules 'see' the fully 
atomistically resolved molecules as coarse-grained 
molecules. Hence, their intermolecular interaction 
is defined by the effective pair potential U cm . To 
ensure that the center-of-mass dynamics governed 
by Eqs. JSJ) and @ of a given explicit or hybrid 
molecule a is correct, the total intermolecular force 
F"g between the atomistically resolved molecule 
a and a coarse-grained molecule /3 is distributed 
among the explicit atoms of the molecule a as 

Fia/3 = ^ — — — F^o, (13) 

where Fi Q/ 3 is the force imposed on the ex- 
plicit tetrahedral atom ia by the coarse-grained 
molecule [3. The explicit tetrahedral atoms in a 
given hybrid molecule interact with other explicit 
atoms in neighboring explicit and hybrid molecules 
through atomistic forces, while the then massless 
center-of-mass particle serves as an effective poten- 
tial interaction site. The total force on a hybrid 
molecule is then according to Eq. J7J a normalized 
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linear combination of atomistic and effective pair 
forces. Using 



T^atom 

iajff ~ 

Fcm 
a/3 ~ 

we obtain from Eq. J7J) 



j/3ia 

Fcm 
Pa 



F a/3 - F 



/3a- 



(14) 
(15) 



(16) 



The force definition in Eq. Q) hence satisfies New- 
ton's Third Law. 

Recall that the effective intermolecular potential 
is determined in such a way that the equations of 
state for the ex and eg models match around the 
state point considered. Therefore, following the 
scheme, as given by Eq. Q, Eq. J3J is implicitly 
satisfied, due to which spurious fluxes are avoided 
at the boundary between the atomistic and coarse- 
grained regimes. 

To summarize, the new adaptive resolution 
scheme for the hybrid atomistic/mesoscale MD 
simulations is a two-stage procedure: 

1. Derive the effective pair potential U cm be- 
tween coarse-grained molecules on the basis 
of the reference all-atom system. 

2. Introduce the interface layer containing the 
hybrid molecules that have additional inter- 
action sites positioned at their centers of 
mass. Define a weighting function w by Eq. 
(Illfl and use Eq. (J7J) for the definition of the 
intermolecular pair forces. Allow molecules 
to adapt their level of resolution according 
to their position in the system as explained 
in the second paragraph of this subsection. 

Finally, since the switching of the resolution can 
be considered as a first order phase transition, 
the adaptive resolution scheme must necessarily 
be employed in combination with a thermostat. 
Because the latent heat is generated in the tran- 
sition regime it is important to use a thermostat, 
which couples locally to the particle motion, e.g., 
Langevin or Dissipative Particle Dynamics (DPD) 
thermostats I4CJ. 



III. COMPUTATIONAL DETAILS 

Temperature Calculation 

In order to treat all n molecules equally regard- 
less of the level of detail, we define, using the 



equipartition principle, the temperature of a sys- 
tem only from translational degrees of freedom. 
The 'translational' temperature 



T = 



2E crn 

3nkn 



1 



3nki 



M a V 2 a , (17) 



when averaged gives the temperature of the system 
|4l) . Here E^" 1 is the total translational kinetic en- 
ergy of the system. We also checked the "particle" 
temperature in the explicit atomistic region. As to 
be expected we find the same temperature. 



Pressure Calculation 

For the same reason also the pressure calcula- 
tion is based on the molecular instead of atomic 
interactions: 



1 

V 



nksT 



3 ^ 

a<(3 



R-a/3 ' F Q /3 



(18) 



where F a p is given by Eqs. l(7)l 142^ . Moreover, 
using Eq. (|18|) for pressure evaluation has two 
additional advantages compared to the pressure 
calculation based on atomic interactions: first, 
Eq. I|18(l is also valid in the case that forces on 
atoms involve internal non-pairwise-additive con- 
tributions, and second, even if all interactions are 
pairwise additive (as in our case), the pressure cal- 
culation based on atomic interactions introduces 
additional fluctuations in the pressure 01 • 



Multiscale Simulation Details 

For computational convenience during the 
proof-of-principle stage we replace our ex-cg model 
with a model in which the whole box contains ex- 
clusively hybrid molecules with four explicit atoms 
and a center-of-mass interaction site. This has 
some technical advantages during the tests of the 
method. The true level of detail of the molecules 
is then determined from the value of w(X a ). For 
later large scale, however, of course the original ex- 
cg model composed of the explicit, coarse-grained, 
and hybrid molecules will be employed. Using 
Eq. 1|11[> in the replacement model for distinguish- 
ing between explicit, coarse-grained, and hybrid 
molecules of our original ex-cg model we can cap- 
ture all definitions of pair forces for all classes of 
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particles in our system in a single expression as 



" iaj0 

w(X a )w(Xf,)Fg$ + [1 - w{X a )w{X p )]5 a ^¥f a ^ 



+ [1 - w{X a )w{Xp)} 



m la m jf 3 



Eia TO » Ej/3 TO J/3 



Fcm 
a/3. 



(19) 



where Fiaj/3 is the total pair force between the ex- 
plicit atom ia of the molecule a and the explicit 
atom j/3 of the molecule (3 and 5 a ^ is the Kro- 
necker symbol. Summing for a ^ f3 



F a/3 = 



(20) 



we regain the total force between molecules a and 
/3 given in Eq. (7J. From Eqs. JSJ, and 
(|19|l follows that in the replacement model with 
only hybrid molecules a hybrid molecule experi- 
ences only translational kicks from other molecules 
in the coarse-grained region (w = 0) and hence its 
center of mass moves exactly as the respective one- 
particle coarse-grained molecule in the original ex- 
cg model. Similarly, in the explicit region (w — 
1) a hybrid molecule experiences only atomistic 
forces and hence its explicit atoms move exactly 
as the explicit atoms in the respective tetrahe- 
dral molecule in the original model. Therefore, the 
model containing only hybrid molecules interact- 
ing via the pair force defined by Eq. 11911 together 
with the applied Langevin thermostat [15j acting 
on each particle in the system (to assure that 
the atom velocities are thermalized in accordance 
with the equipartition principle) exactly mimics 
the original ex- eg model in which the temperature 
would also be held constant by the Langevin ther- 
mostat. From the methodology development point 
of view, these two models are therefore identical. 
This yields the Langevin equation of motion 



d 2 n 

1 dt 2 



m.F^+W^), (21) 
at 



where rrn is the mass of particle i, F.; is the total 
force acting on the respective particle equal to the 
sum of pair interactions given by Eq. (|19f) . T is the 
friction constant, and W is the random force of a 
heat bath^3|. We sample the random force from a 
uniform distribution, since it has been shown that 
there is no advantage of using Gaussian noise for 
the Langevin thermostat [i^. 

The value of the friction constant used in 
our simulations is T — 0.5t _1 where r = 
(e/moo- 2 ) -1 / 2 . The equations of motion are in- 
tegrated for each particle of the system using 



the velocity Verlet algorithm with a 0.005r time 
step. Here, again only for the purpose of test- 
ing the method, we use only one time step in 
the whole system. In the coarse-grained regime 
actually a significantly larger time step could be 
used. Therefore ultimately one would like to in- 
troduce a multiple time step algorithm. Simu- 
lations are performed at temperature T = e/ks 
and number density p = n/V = 0.1/cr 3 . Here 
n = 5001 is the number of molecules in the sys- 
tem, which can be either the explicit, coarse- 
grained or hybrid. If we roughly estimate the 
excluded volume diameter of the coarse-grained 
molecule acG as the distance, where the repul- 
sive effective pair potential between the coarse- 
grained molecules in our simulations equals fc^T, 
i.e., ocg ~ 1-7<t, then the number density p = 



0.1/cr 3 



0.1(a C G/vf/<T C G 



Q.5/a 3 CG 



sponds to a medium dense liquid, which is due 
to the soft effective repulsive interactions rather 
weakly correlated |4^. Periodic boundary condi- 
tions and the minimum image convention |4l| are 
employed. The interaction range in the system is 
given by the range of the effective pair potential 
between molecules and the geometry of tetrahe- 
dral molecules, i.e., the most outer atoms of two 
tetrahedral molecules with centers of mass slightly 
less than 2.31<7 apart still experience the effective 
potential contribution in Eq. i|19|) . Hence, the 
actual interaction range in the system is approxi- 
mately 3.5cr. 

All molecules are initially randomly placed in 
a cubic box of size a = 36.845c. To remove the 
overlaps between them, a 50t long warm-up run 
is performed during which the repulsive interpar- 
ticle potential is capped (see capped interactions 
in Ref. ^J)- Thus, at all interparticle distances, 
which would lead to larger forces between parti- 
cles than a prescribed maximal force, the forces 
defined by the original repulsive pair potential are 
replaced by repulsive central forces of the maximal 
force magnitude. The latter is gradually increased 
from 20e/<7 to 110e/<7 during this warm-up phase. 
Afterwards an additional 250r equilibration run is 
carried out where we set the maximal force magni- 
tude to 10 9 e/cr, which corresponds to interparticle 
distance of 0.27o\ The chosen maximal force mag- 
nitude value is so high that it has no effect on the 
dynamics of molecules in the atomistic region be- 
cause there atoms never come this close together. 
Therefore, by this force capping we only prevent 
possible force singularities that could emerge due 
to overlaps with the neighboring molecules when 
a given molecule enters the interface layer from 
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the coarse-grained side as explained in the previ- 
ous section. Production runs with the 10 9 e/<7 force 
capping are then performed for 7500t, storing con- 
figurations at each 5r time interval for analysis. 
We performed all our MD simulations using the 
ESPResSo package 46], developed at our institute. 

The following reduced Lennard- Jones units[4l| 
are used throughout: m* = m/rriQ, r* = r/cr, 
V* = V/a 3 , T* = k B T/e, U* = U/e, p* = pa 3 /e, 
p* = n/V*, t* = t/r, D* = Dy/rmje/a, where 
D is the self-diffusion constant. Note that in our 
simulations all atoms have a mass m* = 1 while 
all molecules have a mass M* = 4. 



with parameters 7* = 7/e = 0.105, k* = k<t = 2.4, 
and 7-q = tq/ct = 2.31. As one can see from 
figure |3 the obtained effective potential is softer 
than the underlying repulsive interatomic interac- 
tion potential given by Eqs. QJ and J2J since it 
varies more slowly with the interparticle distance. 
This is a general feature of effective potentials for 
polyatomic molecular systems |33j. 

The obtained RDF cm s of the ex and eg systems 
at the temperature T* = 1 and number density 
0.025 < p* < 0.175 using for all number den- 
sities the same effective potential given by Eq. 
<|22|) for the pair interactions between the coarse- 
grained molecules are depicted in figure 0] The 



IV. RESULTS AND DISCUSSION 

Determination of the Effective Potential 

We have determined the effective nonbonded 
pair potential Jj cm * between coarse-grained 
molecules illustrated in figure 03 by using the po- 
tential of mean force (PMF ex ) of the ex system, 
Eq. J3J), at very low number density p* — 0.0025 
as the initial guess. Then we further adjusted it to 
obtain the adequate agreement between RDF cm s 
of the ex and eg systems at the p* =0.1. We have 




FIG. 3: The effective pair potential U cm * , Eq. J23, 
between the coarse-grained molecules, where the po- 
tential of mean force PMF e i of the explicit system at 
p* = 0.0025 and T* = 1 was used as the initial guess. 
The presented function U cm * was determined in such 
a way that the RDF cm s of the explicit (ex) and coarse- 
grained (eg) systems match at the p* — 0.1. 
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VH 7 *{l-exp[- K *(r*-r*)]} 2 (22) 



FIG. 4: Center-of-mass radial distribution functions 
of the explicit (ex) and coarse-grained (eg) systems at 
the temperature T* = 1 and number density 0.025 < 
p* < 2.25. 

RDF cm s are calculated in the range r* £ [0, 5] 
with Ar* = 0.05. From results presented in fig- 
ure 01 it can be observed that although the effec- 
tive potential given by Eq. Q22[) was parameterized 
at the number density p* = 0.1 the agreement of 
RDF cm s of the eg system with the corresponding 
reference RDF cm s of the ex system is good also for 
the lower number densities. In fact, due to weaker 
many-body interactions the agreement is even bet- 
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ter for the systems with lower density. Since the 
effective potentials are density dependent [U, to 
obtain a better agreement for higher densities the 
effective potential should be reparameterizedj^. 
Probably the functional form of the effective po- 
tential will also change due to increased contribu- 
tion from many-body interactions. 

As a quantitative measure of accuracy of the 
evaluated RDF cm s we define a penalty function f p 
as 



g™(r*)) 2 eM-r*)dr*, (23) 



where <7g™, which is taken as a reference, is the 
RDF cm of the ex system. The function exp(— r*) 
is employed to pen alize more strongly deviations at 
small distances 34| . The values of f p for eg systems 
with the number density 0.025 < p* < 0.175 are 
reported in table |U As expected the f p grows with 



p* 


f P ■ w 6 


0.025 


0.0401 


0.075 


0.1101 


0.1 


0.2067 


0.125 


0.4361 


0.150 


0.9257 


0.175 


1.9025 



TABLE I: Penalty function f p defined by Eq. as 
a function of number density p* for RDF cm s gc™(r*) 
of the coarse-grained systems in which particles are 
interacting via the effective potential £/ cm * given by 
Eq. The RDF cm s g c e ™(r*) of all-atom systems 

at the corresponding p* are taken for the reference 
RDF cm s. 

the growing density of the system. 

From the RDF cm we can evaluate the average 
number of neighbors of a given molecule within a 
sphere with the radius r* as 

n n (r*) = p*[ g(r*)Airr* 2 dr*. (24) 
J o 

The n n (r*) for the ex and eg systems at the tem- 
perature T* = 1 and number density 0.025 < p* < 
0.175 are shown in figure [SJ Despite the devia- 
tions between the RDF cm s of the ex and eg sys- 
tems the corresponding average numbers of neigh- 
bors exactly match, indicating, together with the 
RDF cm s presented in figure that the eg model 
with the molecules interacting via the effective po- 
tential given by Eq. H22fl reproduces well the struc- 
ture of the underlying ex system. 

To show that the eg system with the effective 
potential from Eq. (|22|l is at the same state point 
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FIG. 5: The average number of neighbors n„(r*) of 
a given molecule as a function of distance for explicit 
(ex) and coarse-grained (eg) systems at the tempera- 
ture T* — 1 and number density 0.025 < p* < 2.25. 



as the original ex system at the same tempera- 
ture and density we also evaluated the pressure in 
the system. The equations of state for the ex and 
eg models are shown in figure EJ The resulting 




0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 
P 

FIG. 6: Pressure p* in the explicit (ex) and coarse- 
grained (eg) systems at the temperature T* — 1 as a 
function of the number density p* of the system. 
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equations of state are similar to the case of col- 
loidal sphere systems [3^. This similarity is to be 
expected since the tetrahedral molecule as well as 
the spherical coarse-grained molecule can be con- 
sidered as spherical colloidal particles with a hard 
core and a soft coating layer. From figure we can 
conclude that the eg model with the particles in- 
teracting via the effective potential from Eq. (|22ll 
reproduces the equation of state of the underly- 
ing ex system up to p* = 0.125, at which the two 
pressure curves start deviating from each other, 
reflecting the fact that the effective potential is 
density dependent. In order to minimize the arti- 
facts originating from our parameterization of the 
effective potential while still simulating the liquid 
state we have decided to perform all our MD simu- 
lations of the hybrid atomistic/mesoscopic (ex-cg) 
model at the state point with T* = 1 and density 
p* = 0.1. 



Statistical Properties 

The new adaptive resolution scheme is tested by 
comparing the computed statistical properties of 
the ex-cg model with the corresponding properties 
of the reference fully atomistic ex system. 

Figure d (a) displays the RDF cm s calculated 
from center-of-mass positions of all molecules in 
the box of the ex and ex-cg systems at p* = 0.1 
and T* = 1. Shown is the case with the width of 
the interface layer 2d* — 2.5. Depicted are also 
the corresponding local RDF cm s for the explicit 
{ex-cg / ex) and coarse-grained regions {ex-cg / eg) 
of the ex-cg model. As in all simple liquids, e.g., 
methane, the RDF cm s are zero at short distances 
between molecules' centers of mass because repul- 
sive forces prevent overlapping of molecules. Then 
the functions increase rapidly to the first peak. 
With increasing distance RDF cm s reach the lim- 
iting value 1 after few oscillations, indicating that 
there is no order at long distances. The average 
numbers of neighbors n n {r*) of a given molecule as 
functions of distance are illustrated in figured (b). 
The number of nearest neighbors in the first layer 
corresponding to the first minimum in the RDF cm 
is about 11. For comparison, the corresponding 
experimental value from X-ray diffraction for liq- 
uid methane at T = 92 K is approximately 12 47] . 

All calculated RDF cm s of the ex-cg model and 
n n {r*) are in good agreement with the reference 
RDF cm and n n {r*) of the ex model indicating that 
the structure of the underlying all-atom system 
is well reproduced using the adaptive resolution 
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FIG. 7: (a) Center-of-mass radial distribution func- 
tions for all molecules in the box of the all-atom {ex) 
and hybrid atomistic/mesoscopic {ex-cg) systems at 
p* — 0.1 and T* — 1. Shown are also the corre- 
sponding center-of-mass radial distribution functions 
for only the explicit molecules from the explicit region 
{ex-cg/ex) and for only the coarse-grained molecules 
from the coarse-grained region {ex-cg/cg). The width 
of the interface layer is 2d* = 2.5. (b) The correspond- 
ing average numbers of neighbors n n {r*) of a given 
molecule as functions of distance. The different curves 
are almost indistinguishable. 



scheme. This is further confirmed by the values of 
the penalty function f p reported in table [HJ evalu- 
ated for different interface layer widths. From the 
results given in table ITT1 we can see that the local 
structure in the explicit region of the ex-cg model 
is exactly reproduced independently of the inter- 
face layer width. In contrast, the accuracy of the 
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2d* 


f P • io J 


ex-cg 


ex-cg/ex 


ex-cg/cg 


2.5 


0.0821 


0.0009 


0.2109 


3.0 


0.0893 


0.0008 


0.2158 


4.0 


0.1058 


0.0016 


0.2197 


5.0 


0.1324 


0.0023 


0.2315 


6.0 


0.1583 


0.0033 


0.2367 


7.0 


0.1944 


0.0047 


0.2449 


8.0 


0.2301 


0.0076 


0.2616 


9.0 


0.2695 


0.0082 


0.2489 


10.0 


0.3159 


0.0101 


0.2632 



TABLE II: Penalty function f p defined by Eq. as 
a function of the interface layer width 2d* for RDF cm s 
g C eT-c 9 (r*), 9e?- cg/ex {r*), g c ™- cg/cg {r*) of the hybrid 
atomistic/mesoscopic model at p* — 0.1 and T* = 
1- geT-c g ( r *) is the RDF cm of all molecules in the 
box where all molecules are considered indistinguish- 
able, gl™- cg / ex (r*) is the RDF cm of only the explicit 
molecules from the explicit region while g^- cg / cg (r*) 
is the RDF cm of only the coarse-grained molecules 
from the coarse-grained region. The RDF cm ge™{r*) 
of all-atom system at the corresponding p* and T* is 
taken for the reference RDF cm . 



local structure reproduction in the coarse-grained 
region of the ex-cg model depends on the accuracy 
of the effective potential parameterization in the 
eg model (see the value of f p given in table [I] for 
p* = 0.1). The computed f p s for the total RDF crra 
of the ex-cg model lie up to 2d* — 8.0 in between 
the values for the local explicit and coarse-grained 
RDF cm s. 

For narrow interface layers with 2d* < 2.5 the 
adaptive resolution scheme fails to work. The 
molecules are not given enough space and time to 
adapt their orientation to their environment and 
the system is not properly equilibrated in this case. 
The value 2d* — 2.5, for which the adaptive res- 
olution scheme gives the acceptable results, can 
be rationalized by the fact that the interface layer 
width should at least exceed the maximal range 
of interaction, which is the range of the effective 
potential, namely 2.31er. 

To demonstrate that the thermodynamic prop- 
erties are correctly reproduced by the adaptive res- 
olution scheme the temperature and pressure of 
the system as a functions of the interface layer 
width are given in tables ITTT1 and llVl respectively. 
The calculated temperature profile reported in ta- 
ble [ffi] shows that the system is at the right tem- 
perature and that all degrees of freedom are prop- 
erly equilibrated in accordance with the equipar- 
tition principle regardless of the interface layer 
width. The results given in table IIVI also show 



that the adaptive resolution scheme succeeds in 
reproducing the pressure of the underlying fully 
atomistic system. 

In order to check that the chemical potentials in 
the atomistic and coarse-grained regions are equal 
as required by the condition we report in ta- 
ble [V] the average number of molecules in different 
regions of the system. In table E] we also give the 
number of degrees of freedom nno f in the system 
defined as 



n D OF = 3 } [w a N + (1 - w a )}, 



(25) 



where TV = 4 is the number of explicit tetrahe- 
dral atoms in a molecule and w a is the value of 
the weighting function defined by Eq. I|ll|) for the 
molecule a. The summation in Eq. 125|) goes over 
all n molecules of the system. The results show 
that there is no molecule number bias in the sys- 
tem and that the ex-cg system is homogeneous. 
Note that ti^of is greatly reduced employing the 
adaptive resolution scheme on the ex-cg model in 
comparison with the fully atomistic model. The 
time evolution of the number of molecules in dif- 
ferent regions of the system together with the 
nDOF is illustrated in figure [SJ The results clearly 
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FIG. 8: Time evolution of number of molecules in a 
explicit (n e x) , coarse-grained (n C9 ), and interface re- 
gions (ni n t) in the hybrid atomistic/mesoscopic model 
with the 2.5a interface layer width. Time evolution of 
number of degrees of freedom in the system (udof) is 
depicted in the inset. 

demonstrate that the system is in thermodynami- 
cal equilibrium, indicating that the conditions Q 
are satisfied by our adaptive resolution scheme. 

Although we have parameterized the effective 
potential to reproduce the structural properties of 
the ex system we can also compare the dynamical 
properties of the ex-cg with the ex model. For 
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TABLE III: Average temperature as a function of the interface layer width 2d* . T, Tex, T cg , and Tint are the 
average temperatures of the total system, the explicit, coarse-grained, and interface layer regions, respectively, 
calculated by Eq. 1171 . T%x and T°^l are the average temperatures of the explicit and interface layer regions, 
respectively, calculated from total velocities (translational±vibrational±rotational) of explicit atoms in molecules. 
ex and C9 denote the all-atom and coarse-grained systems, respectively. 
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TABLE IV: Average pressure calculated using Eq. 
118H as a function of the interface layer width 2d* . ex 
and C9 denote the all-atom and coarse-grained sys- 
tems, respectively. 

that purpose we have computed the self-diffusion 
coefficient, which is evaluated from the center-of- 
mass displacements using the Einstein relation 

D* = - lixn (|R " (r) - R " ( ° )|2) , (26) 

where R*(i*) is the center-of-mass position of 
the molecule a at time t* and averaging is per- 
formed over all molecules and all choices of time 
origin. The self-diffusion coefficient of the ex and 
eg models calculated in the microcanonical ensem- 
ble, where the Langevin thermostat is switched off 
after the initial warm-up run, are 0.24 and 0.30, 
respectively. The corresponding values of the self- 
diffusion coefficient for the ex, eg, and ex- eg mod- 
els with the Langevin thermostat switched on are 
0.12, 0.14, and 0.13, respectively. Since we use the 
same time and length scales in all our three mod- 



TABLE V: Average number of molecules as a func- 
tion of the interface layer width 2d*. n ex , n cg , and 
riint are the average number of molecules in the ex- 
plicit, coarse-grained, and interface layer regions, re- 
spectively, udof is the average number of degrees 
of freedom defined by Eq. 1251 . For orientation: in 
the system with 2500 coarse-grained molecules, 2500 
four atom explicit molecules, and no hybrid molecules 
riDOF = 37500. ex and C9 denote the all-atom and 
coarse-grained systems, respectively. 



els different values of the self-diffusion coefficient 
in the ex and eg models (with no Langevin ther- 
mostat applied) indicate that the coarse-grained 
molecules experience a slightly smaller intermolec- 
ular frictional hindrance in their motion compared 
to the explicit molecules. This indicates that the 
effective pair potential given in Eq. (1221 intro- 
duces an effective time scale shift in the coarse- 
grained regime. This is known from other studies 
(TT3- E3| - where one actually takes advantage of that 
in order to reach very long simulation times [4e|. 
The apparent self-diffusion coefficient values of the 
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ex and eg models are lower and much closer to- 
gether when the Langevin thermostat is applied 
due to the frictional forces arising from the cou- 
pling to the thermostat unlike to the case of poly- 
mers. There typically the friction of the thermo- 
stat is negligible compared to the friction between 
monomers. Since the self-diffusion coefficient of 
the ex-cg model is close to the corresponding val- 
ues for the ex and eg models we can conclude that 
the center-of-mass dynamics of the molecules is 
similar in all three models. 

As the final test to demonstrate the reliability of 
the adaptive resolution scheme we have computed 
the number density profile of the ex-cg model. 
The results for the system with 2d* = 2.5 and 
2d* = 10.0 are presented in figures ED (a) and (b), 
respectively. The results in figure[!|](a) for the case 
of 2d* — 2.5 show that the explicit and coarse- 
grained regions have the same homogeneous den- 
sity as the reference system while the density in the 
transition regime undergoes an oscillation around 
the reference value Pq =0.1 with a magnitude of 
approximately 0.05p*,- ln the case of 2c?* = 10.0 
(figure (b)) a 5% drop in the density occurs in 
the transition regime, which is compensated by the 
slight increase of the density in the explicit and 
coarse-grained regions. 

This artifact can be explained by considering 
the results displayed in figure ^| (a) , where the 
pressure of the system containing only hybrid 
molecules as a function of the constant value of 
the weighting function w (corresponding to the sit- 
uation in the interface layer) is illustrated. The 
pressure is increased in comparison to the pres- 
sure in the reference system. Clearly the increase 
is most prominent for the most 'artificial' case with 
w = 1/2, indicating that there still is a small 'pres- 
sure barrier' in the interface region causing the 
density dip. This is also evident from results in 
figure ITUI (b), where the RDF cm and the poten- 
tial of mean force of the system with w — 1/2 are 
shown. The effective potential in a system con- 
taining only hybrid molecules with constant value 
of the weighting function changes in comparison 
to the all-atom system. This means that the hy- 
brid molecules in the interface layers of the ex-cg 
model experience too strong effective interaction 
leading to the pressure variations in the interface 
layer. Since the artifact occurs at constant values 
of w it is an artifact of the linear combination of 
forces in Eq. Q and not of the functional form of 
the weighting function w. It must be emphasized, 
however, that this artifact of the proposed adap- 
tive resolution scheme is within a 5% error and 
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FIG. 9: (a) Normalized density profile in the x direc- 
tion of the hybrid atomistic/mesoscopic model with 
the 2.5a interface layer width. Vertical lines denote 
boundaries between atomistic, coarse-grained and in- 
terface regions of the system, (b) The same as in (a) 
but for the 10. Oct interface layer width. 



that similar artifacts, occurring at the boundary 
of two domains with different level of detail, are 
also characteristic for other hybrid schemes [49|. 

The pressure variations in the interface layer 
could cause a spurious reflection of molecules from 
the boundary. However, the results presented in 
figure 1111 where the time evolution of two diffu- 
sion profiles is monitored for molecules that are ini- 
tially localized at the two slabs with a*/10 width 
neighboring the interface layer, show that this is 
not the case. The molecules initially localized at 
the two slabs spread out symmetrically with time. 
This is because the self-diffusion coefficients of all 
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FIG. 10: Artifacts of the adaptive resolution scheme, 
(a) Average pressure p* in the system containing only 
hybrid molecules as a function of the constant value 
of the weighting function w. (b) Center-of-mass radial 
distribution functions for the explicit system (ex) and 
the system containing only hybrid molecules with w — 
const = 1/2 (ea;-c5(w=const=l/2)) at p* = 0.0025 
and T* — 1. The inset also shows the corresponding 
PMF ex and PMF el _ cs(m=cmsfcl / 2 ) determined from 
the systems with p* = 0.0025 using Eq. 



models are approximately the same in the case of 
the applied Langevin thermostat. Thus, the two 
distributions occupy at time t* regions with mean 
square radius 



(\x*{t*) -x*(0)| 2 ) ~ 2D*t* 



(27) 



where x*(0) is the center of the distribution at 
time t* = 0, which is equal to — d* — a*/20 for 
the left slab (figure El (a)) and d* + a*/20 for 



FIG. 11: Time evolution of diffusion profiles for the 
molecules that are initially, at time t* = 0, localized 
at two neighboring slabs of the mid interface layer with 
2d* =2.5 (n is the number of this molecules with the 
center-of-mass position at a given coordinate x*). The 
width of the two slabs is a*/10. Vertical lines denote 
boundaries of the interface layer, (a) The diffusion 
profile, averaged over 500 different time origins, at t* = 
0, t* — 10, and t* = 50 for the molecules that are 
initially localized at the slab on the coarse-grained side 
of the interface region, (b) The same as in (a) but for 
the molecules that are initially localized at the slab on 
the atomistic side of the interface region. 



the right slab (figure ITT1 (b)). Since the diffusion 
profiles are symmetrical at any given time we can 
conclude that the artifact at the interface layer is 
too small to have any severe effect on the diffusion 
of molecules across the interface layer. 
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V. CONCLUSIONS 

A novel approach for efficient hybrid atom- 
istic/mesoscale molecular dynamics (MD) simula- 
tions has been presented in this paper. The new 
adaptive resolution MD simulation scheme dynam- 
ically couples the atomic and mesoscale length 
scales of the studied system by allowing an on- 
thc-fly dynamical interchange between molecules' 
atomic and mesoscopic levels of description. In our 
approach the number of degrees of freedom is al- 
lowed to fluctuate during the course of simulation 
while the statistical properties of the underlying 
all-atom model are properly reproduced for both 
levels of detail. Since the purpose of this paper was 
to develop the method, we restricted ourselves to 
a most simple model system of a molecular liquid 
with short-range repulsive interactions. Even on 
this level a number of questions still have to be 
tackled, as there is the application to higher liq- 
uid densities for complex molecular liquids or the 
proper analysis and understanding of the the time 
scale problem, when the diffusion on the coarse- 
grained level is faster than on the all-atom level. 
For separate runs at the two levels of detail this is 
understood, while the occurrence within one sim- 
ulation box still poses some conceptual problems. 
On the other hand the current density is in be- 
tween a typical small molecule liquid and a poly- 
meric fluid. Thus we expect that already this ap- 
proach can be generalized and applied to differ- 



ent realistic soft condensed matter systems where 
both atomic and mesoscopic length scales have to 
be considered. This can be either polymer solu- 
tions and melts, such as a synthetic or biological 
macromolecule embedded in a solvent. Similarly 
our method should also find application for other 
polymer systems (same force field), molecular liq- 
uids such as methane (same geometry) or water 
(tetrahedral clusters), etc., enabling to reach much 
larger length and time scales than for all-atom MD 
simulations. In all cases the aim is to treat in a 
simulation only as many degrees of freedom as ab- 
solutely necessary for the question considered. In 
this sense the region of higher detail can be ei- 
ther given by a geometrical constraint, e.g., close 
to a surface, or even be chosen on demand due to 
specific local conformations of a (macro-) molecular 
system. Work along these lines is underway. 
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